Selecting putative drought-tolerance markers in two contrasting soybeans

Identifying high-yield genotypes under low water availability is essential for soybean climate-smart breeding. However, a major bottleneck lies in phenotyping, particularly in selecting cost-efficient markers associated with stress tolerance and yield stabilization. Here, we conducted in-depth phenotyping experiments in two soybean genotypes with contrasting drought tolerance, MUNASQA (tolerant) and TJ2049 (susceptible), to better understand soybean stress physiology and identify/statistically validate drought-tolerance and yield-stabilization traits as potential breeding markers. Firstly, at the critical reproductive stage (R5), the molecular differences between the genotype’s responses to mild water deficit were explored through massive analysis of cDNA ends (MACE)-transcriptomic and gene ontology. MUNASQA transcriptional profile, compared to TJ2049, revealed significant differences when responding to drought. Next, both genotypes were phenotyped under mild water deficit, imposed in vegetative (V3) and R5 stages, by evaluating 22 stress-response, growth, and water-use markers, which were subsequently correlated between phenological stages and with yield. Several markers showed high consistency, independent of the phenological stage, demonstrating the effectiveness of the phenotyping methodology and its possible use for early selection. Finally, these markers were classified and selected according to their cost-feasibility, statistical weight, and correlation with yield. Here, pubescence, stomatal density, and canopy temperature depression emerged as promising breeding markers for the early selection of drought-tolerant soybeans.

. Transcriptomic analysis of MUNASQA and TJ2049 genotypes under drought. Heat-map of all DEGs for MUNASQA and TJ2049 in drought conditions. Scale color indicates green for up-regulation and red for downregulation (a). GO enrichment in MUNASQA and TJ2049 comprises biological processes (BP, in red), molecular function (MF, in blue), and cellular component (CC in green). Relevant categories showing enrichment of DEGs for both genotypes are depicted. GO terms were plotted after applying an FDR = 0.1 Bubble size correlates with enrichment factor values; for each bubble size, the P-value is indicated (b). Venn diagram for all DEGs in MUNASQA and TJ2049 under drought conditions. DEGs were plotted after applying an FDR = 0.1 (c). Validation by qRT-PCR of ten genes selected from RNA-Seq. Log 2 fold change (log2FC) was calculated based on the comparison of drought vs control for each genotype (d). Three biological replicates were used, and the experiment was performed twice with similar results. www.nature.com/scientificreports/ water deficit, and independently of the phenological stage (Table 1). MUNASQA reached maximum activity for all enzymes after 4 d of water deficit, while the highest activity for TJ2049 occurred to 8 d after the stress onset, except for CAT. Noticeably, under non-stressed conditions, all enzymatic activities, excluding CAT, showed higher activity in the tolerant genotype than in the susceptible one, strengthening the hypothesis that TJ2049, compared to MUNASQA, presents delayed stress perception and response mechanisms. In fact, and according to Carvalho et al. 16 , a successful response to water deficit may depend not only on which enzymes are activated but also on the activation timing. Moreover, the general gene expression profile regarding these enzymes regulation was consistent with the biochemical results.
Regarding PRO, one of the most common osmoprotectant in plants 18 , MUNASQA and TJ2049 accumulated the osmolyte in response to water deficit over time and in both V3 and R5 stages (Table 1). Several studies have demonstrated a direct correlation between high osmoprotectant accumulation and drought tolerance in many crops 19 . Here, the tolerant genotype MUNASQA exhibited a higher and more rapid accumulation of PRO after 4 d of water deficit at both phenological stages. In agreement with our results, a recent study in soybean reported higher PRO accumulation in the drought-tolerant genotype A5009 RG, compared with the drought-susceptible ADM50048 20 .
When analyzing MDA production, an indicator of lipid peroxidation and stress severity 21 , a higher and more rapid accumulation was detected in TJ2049 plants in response to water deficit, compared to MUNASQA (Table 1).
Drought also affects leaf pigments content 21 . Changes in photosynthetic pigments can alter various lightharvesting processes, while the accumulation of photoprotective compounds plays an essential role in preventing photo-oxidative damage 22 . Here, MUNASQA and TJ2049 showed alterations in pigment content under water deficit ( Table 1). The CHL was significantly reduced over time due to stress in both genotypes, but this reduction was significantly lower in the tolerant one. Similar results were found in drought-tolerant maize that showed lower CHL reductions under stress than susceptible genotypes 23 . Regarding CAR levels, TJ2049 and MUNASQA showed an increased synthesis in response to drought, although the last one exhibited greater and faster accumulation.
According to our results, MUNASQA drought tolerance is strongly related to a rapid stress-sensing and response capacity and an efficient ROS (reactive oxygen species) scavenging system, both essential mechanisms for stress tolerance in numerous species 24 . Growth and yield markers. The ability to produce high seed yield or biomass under limited water access is considered the optimal indicator of drought tolerance in crops [25][26][27] . We evaluated the effect of mild water deficit in MUNASQA and TJ2049 growth and yield by monitoring changes in various markers related to biomass and seed production, including leaf area index (LAI), leaf area ratio (LAR), the net assimilation rate (NAR), relative growth rate (RGR), crop growth rate (CGR), relative yield and DSI (Table 2).
Total leaf surface area and LAI are strongly associated with canopy interception, evapotranspiration, and photosynthesis 28 . Here, independently of the phenological stage or water availability, MUNASQA plants exhibited a higher LAI compared to TJ2049, indicating a larger assimilatory capacity and, as a consequence, photosynthetic potential. In vegetative stages, a higher LAI denotes a more rapid canopy development, favoring greater and faster soil coverage and thus less water loss from direct evaporation. In general, drought suppresses leaf initiation and growth and consequently affects LAI 29 . Therefore, a decrease in LAI is expected in plants exposed to water scarcity. As expected, MUNASQA and TJ2049 plants showed a reduction in LAI in response to drought, more noticeable after 8 d of stress and in the R5 stage.
Drought also alters the LAR: the leaf area development in relation to the total biomass produced 30 . Here, water deficit affected the LAR in both genotypes, but only during the vegetative stage. The highest ratio between plant leaf area and biomass is reached at the beginning of the plant life cycle 31 , which explains the highest LAR values at the first sampling day during the vegetative stage.
In MUNASQA and TJ2049, the relationship between leaf area expansion and the biomass produced over time (NAR) was also reduced due to water deficit at both phenological stages. Changes in photosynthetic efficiency were more significant at V3, where both genotypes exhibited greater NAR values. Noticeably, wellwatered MUNASQA plants showed a significantly higher NAR than TJ2049 ones, while, in response to stress, MUNASQA's NAR increased in contrast to TJ2049 values, which were reduced. Although not as accentuated, a similar pattern was observed in plants exposed to water deficit in R5. Here, LAI, LAR, and NAR results indicated that, in response to drought, MUNASQA plants regulated photosynthates allocation to leaves and the maintenance of photosynthetic efficiency more efficiently than TJ2049 ones.
In addition, the relative growth rate (RGR) or biomass produced over time was significantly reduced by water deficit in the vegetative stage (Table 2). Interestingly, the reduction in RGR was more pronounced in TJ2049 plants. Regarding CGR, significant differences were observed in the V3 stage after 8 d of drought, while in R5 were detected after 4 d.
These results agree with the differences in yield and yield-DSI exhibited by MUNASQA and TJ2049 after water-deficit treatments in V3 and R5 (Fig. 2). According to these findings, when drought was applied in the V3 stage, TJ2049 showed a distinct but not significant yield penalty and a yield-DSI considerably higher than MUNASQA. Moreover, after a mild water deficit in R5, a highly moisture-sensitive phenological stage, TJ2049 exhibited the largest yield loss and a significantly higher yield-DSI than MUNASQA. Results that also agree with the ones reported by Pardo et al. 13 .
Water use markers. Maintaining tissue/cellular water content and/or metabolic activity at low water potentials are physiological strategies to survive drought 32 . Traits like pubescence, leaf thickness, stomatal den- Table 1. Effect of mild water deficit on stress-response enzymatic markers measured in MUNASQA and TJ2049. SOD, APX, POX, and CAT activities, together with PRO, MDA, CHL, and CAR contents, were obtained from plants submitted to water deficit (Ψs = − 0.65 MPa) and well-watered treatments (Ψs = − 0.05 MPa) applied in V3 and R5 stages. Two independent experiments (n = 10 per genotype/treatment) were conducted, assessing parameters at 0, 4, and 8 d after stress (DAS) imposition. Additionally, 50 plants per genotype and the following treatments: 1: Control, 2: V3-Stress, and 3: R5-Stress, were harvested at physiological maturity to obtain relative yield. Average values followed by the same uppercase letter in the column and the same lowercase letter in the row do not differ statistically among them within each phenological stage, according to Tukey's HSD test at 5%. The strength of association between markers evaluated in V3 and R5 stages (n = 240) and between markers and yield (n = 300) was measured by Pearson's correlation analysis adjusted by Bonferroni (P > 0.05 indicated as ns; P < 0.05 indicated as *; P < 0.01 ** and P < 0.001 ***). Correlation coefficients (r 2 ) were classified as "S: Strong" (> ± 0.60) and "W: weak" (below ± 0.59). Significant values are in bold.      www.nature.com/scientificreports/ Table 2. Effect of mild water deficit on growth markers measured in MUNASQA and TJ2049. LAI, LAR, NAR, RGR and CGR were assessed in plants submitted to water deficit (Ψs = − 0.65 MPa) and well-watered treatments (Ψs = − 0.05 MPa) in V3 and R5 stages. Two independent experiments (n = 10 per genotype/treatment) were conducted, assessing parameters at 0, 4, and 8 d after stress (DAS) imposition. Additionally, 50 plants per genotype and the following treatments: 1: Control, 2: V3-Stress and 3: R5-Stress, were harvested at physiological maturity to obtain relative yield. Average values followed by the same uppercase letter in the column and the same lowercase letter in the row do not differ statistically among them within each phenological stage, according to Tukey's HSD test at 5%. The strength of association between markers evaluated in V3 and R5 stages (n = 240) and between markers and yield (n = 300) was measured by Pearson's correlation analysis adjusted by Bonferroni (P > 0.05 indicated as ns; P < 0.05 indicated as *; P < 0.01 ** and P < 0.001 ***). Correlation coefficients (r 2 ) were classified as "S: Strong" (> ± 0.60) and "W: weak" (below ± 0.59). Significant values are in bold. www.nature.com/scientificreports/ sity, closure, slow wilting, canopy temperature, RWC, and WUE are essential to determining drought tolerance in plants. Thus, all these water-use parameters were assessed in MUNASQA and TJ2049 in response to drought. After 21 d of water deficit, the genotypes exhibited drought-induced adaptations in leaf morphology traits 33 , such as stomata, trichrome density, and leaf thickness ( Table 3). As expected, the stomatal density was substantially altered in response to drought. Here, TJ2049 plants exhibited a considerable decrease, while MUNASQA ones showed an 89 and 65% increase on abaxial and adaxial leaf surfaces. Moreover, no changes in pubescence were observed on TJ2049 plants, while in MUNASQA, the trichrome density was increased, especially in the abaxial leaf surface. In response to drought, plants can reduce their stomatal size, density, or aperture and develop higher pubescence 34 , mainly on the abaxial surface 35 . The increase in MUNASQA's stomatal density could indicate the genotype's ability to produce smaller but denser stomata and, therefore, reduce transpiration by a quicker onset of stomatal regulation. This modification, combined with denser trichomes, could enhance the boundary layer resistance, increase the air's moisture outside the stomata and minimize water loss during drought, precluding significant growth penalties in terms of photosynthetic activity.

Genotype and
Regarding leaf thickness, no significant differences were detected under well-irrigated conditions. However, after 21 d of water deficit, MUNASQA leaves were considerably thinner, while TJ2049 ones were thicker. The knowledge about the links between leaf morpho-anatomy and its function under non-stressed/stressed conditions is relatively poor, especially for traits like leaf thickness 36 that is strongly related to transpiration 37 and reported by some authors as a drought-tolerance trait that maintains turgor pressure and enhances photosynthesis 38 . Yet, this feature, only apparent in TJ2049 plants under water deficit, could not be associated with transpiration adjustments, photosynthesis increase, or any other drought-tolerant feature. . Different letters indicate significant differences at P < 0.05 (two-way ANOVA). Error bars represent SE from independent experiments, n = 300 per trial. Table 3. Effect of mild water deficit on leaf morphology of MUNASQA and TJ2049. LT, TD_AB, TD_AD, ST_AB, SD_AD, and stomatal aperture were assessed in plants submitted to water deficit (Ψs = − 0.65 MPa) and well-watered treatments (Ψs = − 0.05 MPa) in R5 stage (except for stomatal aperture applied in V3). For LT, SD_AB, SD_AD, TD_AB and TD_AD, an independent experiment (n = 5 per genotype/treatment) was conducted, assessing parameters at 3, 10 and 21 days after stress (DAS) imposition. Here we showed the data corresponding to 21 DAS (n = 10 measured per sample). For stomatal aperture, three independent experiments (n = 40 stomatal measurements per genotype/treatment) were conducted, and the stomata evaluation was performed 72 hs after stress imposition. Average values followed by the same uppercase letter do not differ statistically according to Tukey's HSD test at 5%. www.nature.com/scientificreports/ Meanwhile, the regulation of stomatal aperture reinforces the drought-tolerant character of MUNASQA. Under non-stressed conditions, the susceptible TJ2049 showed more opened stomata (⁓22% more than MUNASQA). Moreover, after 3 d of water deficit in V3, this difference was increased to almost ⁓50% of stomatal aperture (Table 3). Although stomata represent a small percentage of the leaf lamina, large amounts of water evaporate through them 33 . Thus, the lack of stomatal control in TJ2049 might explain the fast-wilting phenotype in response to air desiccation (Fig. 3). Measuring the water loss of detached leaves is a selection method for drought tolerance 39 . Here, R5 leaves of both genotypes were removed and air-dried. After 48 h, MUNASQA exhibited a greener, healthier and slow-wilting phenotype, previously linked to drought tolerance in studies with soybean cultivars 17 .
MUNASQA water-saving behavior was also confirmed through parameters such as RWC and WUE, strongly regulated under drought ( Table 4). The RWC, a time-specific measurement of the hydric status of a plant, is considered a physiological character recommended for drought-tolerance selection 40 . In V3 and R5 stages, wellirrigated MUNASQA plants showed higher RWC than TJ2049, even with a smaller canopy. As expected, under water deficit, the RWC values were reduced by 14% in MUNASQA and 20% in TJ2049, indirectly confirming the effectiveness of the stress imposed.
Water use efficiency (WUE), referring to the biomass produced per water unit, has been widely used as a breeding target in many rainfed crops, including soybean. Conservative water-use strategies are associated with high leaf WUE 34 . In agreement, and contrarily to TJ2049, V3 and R5 MUNASQA plants showed a gradual increase of WUE in response to water deficit that agrees with the tighter regulation of stomatal movements and the reduced water loss observed in the genotype (Table 4). Moreover, considering the discrete NAR reduction and the maintenance of a ~ 70% RWC under water deficit, we hypothesize that MUNASQA may display a stomatal control based on partial or total/partial closure intervals, therefore reducing transpiration and saving water through a smaller gas exchange (potential photosynthesis) penalty.
The CTD, regarding plant canopy temperature difference with the surrounding air, is considered a surrogate trait for stomatal conductance and a good indication of plant transpiration rate 41 . As expected, in response to drought, MUNASQA plants evidenced lower CTD values, a finding that supports the stomatal aperture results and strongly suggests the genotype water-saving behavior. Plants with higher stomatal conductance transpire more and thus maintain a cooler canopy 42 . Thus, in TJ2049 stressed-plants, the high and positive CTD confirmed a higher stomatal aperture and transpiration rate that agrees with a water-spender behavior. Moreover, TJ2049 also presented higher transpiration rates in unstressed conditions. Finding that could be evidence of a natural and predisposing difference between tolerant and susceptible genotypes.

Markers selection.
Identifying and exploiting phenotyping markers will improve selection strategies for drought tolerance in legumes crops 17 . However, to successfully implement markers in a breeding program, it is imperative to validate their (i) accuracy, (ii) feasibility, and (iii) strength of association with the desired trait. To further understand the marker's contribution to drought tolerance and yield stabilization in MUNASQA and TJ2049, Principal Component Analysis (PCA) was performed.
A first PCA was conducted for all the morphophysiological parameters evaluated together with absolute yield (Fig. 4). Here, data corresponding to MUNASQA stress was separated from the rest of the treatments and genotype in Principal Component (PC) I, showing the biggest dissimilarity (Fig. 4a). Meanwhile, TJ2049 stress data were separated in PC II. However, when comparing the proportion of variance explained by the different PC (Fig. 4b), the first two only explained 53,48%. Moreover, no clear association was observed between the absolute yield and the rest of the parameters, although PRO and LAR exhibited some positive relation.
In PCAs, an increase in the number of comparable variables will reduce the proportion of variance explained by those variables. Therefore, to discriminate which markers better explain the variability between genotypes and treatments, independent PCAs were conducted using subsets of parameters grouped by biological processes in (i) "stress response", (ii) "growth" and (iii) "water use" categories ( Fig. S3). In addition, markers evaluated at V3 and R5 stages were measured by Pearson's correlation to determine their strength of association between the phenological stages (accuracy) and yield stabilization (desired trait) (Tables 1, 2, and 4).
Biochemical parameters such as enzymatic and non-enzymatic ROS scavengers, leaf pigments, and MDA have been confirmed as adaptive responses to desiccation stress, frequently used for selecting plant genotypes under drought 43 . The PCA results showed clear discrimination of MUNASQA and TJ2049 drought responses (Fig. S3a). Here, the first two principal components (PC) explained 96.6% of total variation (PC1 = 75.0% and PC2 = 21.6%). Data were clustered by irrigation treatment in PC1, suggesting that these markers are indicators of phenotypic plasticity. The CHL was associated with well-irrigated plants, while all enzymes, MDA, CAR, and PRO, were related to drought stress.
All these markers showed high accuracy between phenological stages due to significant (P < 0.001), strong (r 2 over 0.88), and positive correlations (Table 1). However, the correlation with yield showed inconsistent outcomes. PRO, MDA, and CAR were significant and positively correlated with yield. Meanwhile, except SOD, all enzymes exhibited significant and positive correlations that were too variable in the strength of association and could not be linked with a specific genotype or treatment. Thus, we considered these "stress-response" markers suitable for discriminating susceptible/resistant responses during early drought-tolerance screenings in soybean. Still, their use as breeding traits is limited due to the high environmental effect. Interestingly, a report 20 found that PRO and CHL were suitable markers for ranking soybean genotypes in response to drought in vegetative stages (5 days after emergence). At the same time, MDA could be useful during R5 as a sensitivity trait.
In legumes, features like LAI, smaller leaf area, leaf area maintenance, and dry matter partitioning have been used to screen for drought tolerance 17 . Here, in the "growth" PCA (Fig. S3b), the first two PC explained 93.0% of total variation (PC1 = 73.0% and PC2 = 20.0%). In PC1, data were clustered by genotype, suggesting that these www.nature.com/scientificreports/ markers explained differences between MUNASQA and TJ2049 growth responses on their intrinsic genetic variability. Moreover, the markers analyzed were only associated with MUNASQA, LAR, and RGR associated with stressed plants. However, only LAI and LAR showed a weak correlation between phenological stage and yield ( Table 2). Water-saving features like denser leaf pubescence, a higher number of stomata, warmer canopies, RWC maintenance, and increased WUE have been associated with drought tolerance in legumes and applied in droughtresistance breeding 17 . Here, "water-use" markers contributed the most to discriminating drought-tolerant and susceptible responses. In the two PCA performed, one for physiological markers and the other for morphological ones, data were clustered by genotype in PC1; thus, these markers are indicators of genetic variability. In CTD, RWC and WUE PCA (Fig. S3c), the first two PC explained of total variation 87.7% (PC1 = 54.6% and PC2 = 33.1%). Markers RWC and WUE (associated with MUNASQA) and CTD (associated with TJ2049) showed significant (P < 0.001), strong (r 2 over 0.87), and positive correlations between phenological stages (Table 4). Moreover, WUE and CTD were significantly associated with yield, showing positive and negative correlations depending on the water treatments applied (Table 4). In the PCA made with morphological "water-use" markers, the first two PC explained 96.6% of total variation (PC1 = 73.8% and PC2 = 22.8%) (Fig. S3d). Here, pubescence and stomata abaxial density was strongly related to MUNASQA. Although the data were insufficient to execute good correlation analysis, these morphological markers have been demonstrated as clear indicators of watersaving strategies in legumes 34 .
A good drought-tolerance marker linked to yield stabilization must also be accurate, cost-effective, if possible, non-destructive, and easily measurable. Hence, after evaluating marker accuracy (amid phenological stages) and assessing which ones better explained the phenotypic variability between genotypes and treatments, a final selection was performed by cost-feasibility (CF) and statistical weight (SW) rankings. Based on the CF values, the degree of complexity and cost to assess each indicator, nine markers were further selected, encompassing the categories 1 (easy and cheap) and 2 (easy and expensive) ( Table 5). Subsequently, after SW re-selection, based on the percentage of variability among genotypes and treatments explained by each parameter, four markers remained with "High" SW in both PC.
Determine which markers are more suitable to assess drought tolerance often represents a challenge and depends on the researcher's criteria, e.g., whether tolerance is based on yield maintenance or intrinsic mechanisms that ensure plant survival at the expense of productivity. During this research, we assessed numerous parameters associated with plant performance and stress responses, aiming to identify a small group of markers related to yield stabilization, stress tolerance, or both, and if possible, non-destructive and easily measurable. Often, some markers are highly accurate (stable during the plant life cycle) but expensive, laborious, and/or time-consuming (low throughput). Thus, without disregarding the importance of accuracy, we also consider markers that are cost-efficient and can be assessed in a high throughput manner as strong candidates for droughttolerance phenotyping.
During this research, the selected markers were (i) stomatal density on the adaxial and (ii) abaxial leaf surface, (iii) trichrome density on the abaxial side, and (iv) CTD. These four traits were chosen as the most efficient phenotyping markers for drought tolerance due to their high accuracy, strong association to water-saving strategies www.nature.com/scientificreports/ Table 4. Effect of mild water deficit on water-use physiological markers measured in MUNASQA and TJ2049. RWC, WUE and CTD were assessed in plants submitted to water deficit (Ψs = − 0.65 MPa) and well-watered treatments (Ψs = − 0.05 MPa) in V3 and R5 stages. Two independent experiments (n = 10 per genotype/treatment) were conducted, assessing parameters at 0, 4, and 8 d after stress (DAS) imposition. Additionally, 50 plants per genotype and the following treatments: 1: Control, 2: V3-Stress and 3: R5-Stress, were harvested at physiological maturity to obtain relative yield. Average values followed by the same uppercase letter in the column and the same lowercase letter in the row do not differ statistically among them within each phenological stage, according to Tukey's HSD test at 5%. The strength of association between markers evaluated in V3 and R5 stages (n = 240) and between markers and yield (n = 300) was measured by Pearson's correlation analysis adjusted by Bonferroni (P > 0.05 indicated as ns; P < 0.05 indicated as *; P < 0.01 ** and P < 0.001 ***). Correlation coefficients (r 2 ) were classified as "S: Strong" (> ± 0.60) and "W: weak" (below ± 0.59). Significant values are in bold. www.nature.com/scientificreports/ under drought, high cost-efficiency (affordable and easily measured), and non-destructive assessment, therefore ideal for high throughput screening using high-resolution imaging. Furthermore, beyond the contributions to our soybean breeding program, we must highlight the practical applications of the in-depth phenotyping results. Overall, these findings (i) confirmed the effectiveness of the methodologies used during the research (e.g., drought imposition, sampling times), (ii) corroborate its successful application in early phenological stages, and (iii) strength the usefulness of MUNASQA and TJ2049 as model genotypes for genetic mapping studies. In this context, our group developed a segregating population of 280 F6 RIL (recombinant inbred lines) families that will be mapped for drought tolerance with SSR (simple sequence repeats or microsatellites) and SNPs (single nucleotide polymorphism) selected based on the four markers chosen in this research (stomatal and trichome densities and CTD). stage (beans beginning to develop at one of the four uppermost nodes with a wholly unrolled leaf) was assessed through transcriptional and leaf morphology analysis. Subsequently, comparative studies were performed to determine the genotype's response to water deficit imposed in V3 and R5 stages. Next, all markers assessed were analyzed according to their strength of association between phenological stages and yield, then were ranked by statistical weight and cost-feasibility (Fig. S1). Irrigation treatments. The volumetric water content (VWC) of each pot was estimated according to Pereyra-Irujo et al. 45 , and the relationship between VWC and water potential (Ψs) was determined 46 . Pots were maintained at 22% of VWC (Ψs = − 0.05 MPa) through daily watering until stress onset. According to Pardo et al. 13 , the water deficit was applied by maintaining the pots at 14% of VWC (Ψs = − 0.65 MPa) for ten days. The desired Ψs was reached in 2-3 days. At the end of stress, plants were fully watered until harvest. The Ψs was daily monitored and recorded. Corrections for soil water status were made by weighing two plants per genotype and treatment every 3 days. The plant water status was monitored through the RWC 47 to ensure stress occurrence.

Plant material and growth conditions.
All drought experiments were carried out for three consecutive years, always applying the previously described irrigation treatments. Water deficit in V3 and R5 phenological stages was imposed in independent plant sets; the sections below detailed the sampling process.

Experiments. MACE-transcriptomic analysis and validation.
Three biological replicates per treatment (Control and Stress) and genotype (MUNASQA and TJ2049) were collected from R5 plants after 72 h of water deficit (n = 12), and RNA from fully mature expanded leaf between nodes 5 to 7 node was isolated for transcriptional analysis.
MACE-Seq libraries and sequencing were performed on an Illumina NextSeq500 machine (1 × 75 bp reads). The conversion was made with bcl2fastq2 software (version 2.19.1), and the cleaning of duplicate sequences was performed with "TrueQuant". In MACE-seq, the TrueQuant barcodes each DNA molecule before PCR amplification. As each barcode-template combination is statistically unique, PCR-duplicates can be identified and eliminated from the dataset to prevent PCR bias. Bases with low sequencing quality were clipped. Next, reads were mapped into genome version "Gmax_275_v2.0.fa" of soybean downloaded with standard parameters from Phytozome.net and Bowtie2 (version 2.2.4). Then, expression analysis was performed by in-house scripts and DESeq2 (R-package).
DEGs were defined at an FDR of 0.1 and listed as either up or down-regulated. A heat map plot was generated using R software (version 3.4.1). Then, hierarchical clustering was applied by considering a cut-off threshold of 8 expression profiles (clusters). Venn diagrams were depicted using the VennDiagram R package.
A GO enrichment analysis was performed using the topGO R package 48 , while the GO annotation file was extracted from agriGO website 49 . Each GO term, containing at least two DEGs, was analyzed by Fisher's exact www.nature.com/scientificreports/ test. The resulting P values were corrected by FDR multiple testing approach. GO terms with an FDR lower than 0.1 were considered for further analysis. Ten DEGs were randomly selected and measured by qRT-PCR assays (Applied Biosystems) to validate MACE results. F-BOX gene was used as an internal reference to standardize the expression of target genes, and the ratio between treatments was calculated according to 50 . All primers used are listed in Suppl. Table 6. Data analysis and primer efficiencies were obtained using LinReg PCR software 51 . Relative expression ratios and statistical analysis were performed using fgStatistics software interface 52 . The cut-off for statistically significant differences was set as P < 0.05, indicated as *.
Leaf morphology measurements. Changes in leaf thickness (LT), adaxial and abaxial stomatal and trichrome densities (SD_AD, SD_AB, TD_AD and TD_AB) were assessed in leaves between nodes 4 to 7 of MUNASQA and TJ2049 R5 plants after 3, 10 and 21 d of water deficit. Five samples per genotype and treatment (n = 20) were taken and fixed in FAA (10% formalin, 5% acetic acid, 50% ethyl alcohol). Diaphanised sections of the central leaflet were used for superficial observations. Different standard colorations were applied according to D' Ambrogio de Argüeso (1986) 59 . Staining samples were visualized in a Leica DM500 optical microscope and photographed with an Arcano (5 Mpx) camera (10 measurements per sample, n = 200).
Stomatal aperture measurements. According to Gudesblat et al. 60 , stomatal apertures were measured in three independent assays, using MUNASQA and TJ2049 plants submitted to 72 h of water deficit in the V3 stage. The aperture of 40 stomata per treatment and genotype (n = 120) was measured in each experiment.
Wilting air desiccation assay. Response to air desiccation was evaluated in the R5 stage. Three whole leaves per genotype (n = 6) were collected and exposed to air desiccation at 32 °C. After 0, 6, 24, 36 and 48 h of air exposure, plants were photographed with a Canon Power Shot SX520 HS (14 Mpx), and the wilting rate was assessed.
Comparative analysis of genotypes responses to water deficit in V3 and R5. The responses of MUNASQA and TJ2049 to water deficit, applied in V3 and R5 stages, were compared by measuring morphophysiological and biochemical parameters grouped by biological processes (BP) ( Table 6). Four treatments were defined: (i) Control-V3, (ii) Control-R5, (iii) Stress-V3 and (iv) Stress-R5, and three sampling times were performed (0, 4 and 8 Table 6. Markers evaluated in MUNASQA and TJ2049, clustered by biological processes (BP). www.nature.com/scientificreports/ d of water deficit). Ten plants per genotype, treatment and time were collected and used for markers evaluation (n = 240). Additionally, for treatments (i), (iii) and (iv), 50 plants per genotype were harvested at physiological maturity (n = 300) to quantify relative yield and calculate the relative yield DSI (Drought Susceptibility Index) according to Fischer and Maurer 61 .
The markers evaluated and their methodologies are detailed below.
As growth indicators, the plant total leaf area (TLA) and biomass (plant total dry weight) were determined. Then leaf area index (LAI) and leaf area ratio (LAR) 66 , the net assimilation rate (NAR) 67 , relative growth rate (RGR) 68 and crop growth rate (CGR) 69 were calculated.
Finally, as water-use parameters, plant RWC and WUE 70 were calculated. Moreover, the canopy temperature was monitored and recorded using a FLIR ONE-3 thermal camera (0.3456 Mpx) to calculate canopy temperature depression (CTD) 71 .
Univariate analysis. Data from stomatal apertures were subjected to a two-way ANOVA (Factor 1: genotype, Factor 2: treatment). The remaining data were analyzed through ANOVA with post hoc contrast by Tukey's HSD test. Data were analyzed with InfoStat statistical package 52 and presented as the arithmetic mean ± SE. Means were considered significantly different at P < 0.05.
All markers were submitted to a PCA to discriminate main associations between markers, genotypes and treatments. However, in PCAs, an increase in the number of comparable variables will reduce the proportion of variance among treatments explained by those variables. Therefore, all markers were grouped by biological processes in (i) "stress response", (ii) "growth", and iii) "water use" sets and subjected to independent PCAs to discriminate which markers better explain the variability between genotypes/treatments. Additionally, the markers were ranked by CF and SW. The markers CF, in terms of their complexity and evaluation cost, was assigned according to 4 categories: easy and cheap (1), easy and expensive (2), complicated and cheap (3) or complicated and expensive (4). Meanwhile, the SW was obtained from PCA variables coefficients (autovectors e1 and e2) that were ranked and classified in "Low" (Low = [− 2, 2]) and High (High = R − [− 2, 2]), according to their weigh on PC1 and PC2. Markers strongly correlated between phenological stages, if possible, with yield, together with CF values of 1 or 2 and "High" SW in both PC, were selected as the most efficient phenotyping markers.